library(mellonMisc)
library(lme4)
library(lfe)
library(mixtools)
library(ineq)


simDiffsMods <- function(n.sims = 100, n.promises = 250, 
                         group.diff = 3, group.size = 5, obs.prom = 200, resp.sd = 1.7) {
  
  sim.output <- dtf(sim = 1:n.sims, fe.sd = NA, re.sd = NA,
                    true.sd = NA, re.clsd = NA, cluster.sd.ratio = NA,
                    re.cl.sd.ratio = NA)
  
  re.cl.wt <- fe.wt <- re.wt<- true.wt <- re.cl.eff <- fe.eff<- re.eff <- 
    true.eff <- matrix(nrow = n.sims, ncol = n.promises)
  
  for(sim in which(is.na(sim.output$fe.sd))) {
    print(sim)
    true.values <- rnorm(n.promises, sd = 0.5, mean = 0)
    true.values[1:group.size] <- true.values[1:group.size] + group.diff
    
    
    true.weights <- prop.table(exp(true.values / sd(true.values)))
    
    
    n.respondents <- (obs.prom * length(true.values)) / 8
    
    respondents <- sample(inverse.rle(list(values = 1:n.respondents,
                                           lengths = rep(8, n.respondents))))
    
    resp.eff <- rnorm(n.respondents, sd = resp.sd)[respondents]
    
    obs.vals <- as.list(rep(NA, length(true.values)))
    
    
    for(ii in 1:length(true.values)) {
      obs.vals[[ii]] <- rnorm(obs.prom, mean = true.values[ii], sd = 3)  
    }
    
    data <- dtf(promise = inverse.rle(list(values = 1:n.promises, lengths = rep(obs.prom, n.promises))), 
                y = unlist(obs.vals), respondents= respondents)
    
    data$y <- resp.eff + data$y
    
    
    prom.res <- lmer(data = data, y ~ (1|promise) + (1|respondents))
    su.re.basic <- summary(prom.res)
    
    re.prom <- ranef(prom.res)$promise[, ]
    
    re.weights <- prop.table(exp(re.prom / attributes(su.re.basic$varcor$promise)$stddev))
    
    data$z <- data$y - mean(data$y)
    
    
    fe.prom <- felm(z ~ 0 | promise + respondents, data  = data)
    
    fe <- getfe(fe.prom)
    fe <- fe[1:length(true.values), "effect"]
    fe.weights <- prop.table(exp(fe / sd(fe)))
    
    
    mix.mod <-  try(normalmixEM(re.prom, maxit = 10000))
    
    
    if(class(mix.mod)!="try-error") {
      cluster.prob <-mix.mod$posterior
      data$cluster <- cluster.prob[data$promise, 1]
      prom.res.cl <- lmer(data = data, y ~ cluster + (1|promise) + (1|respondents))
      
      re.prom.cl <- ranef(prom.res.cl)$promise[, ]
      su.re <- summary(prom.res.cl)
      
      
      re.prom.cl <- re.prom.cl + cluster.prob[, 1] * su.re$coefficients["cluster", "Estimate"]
      
      comb.sd <- sqrt(sd(data$cluster * su.re$coefficients["cluster", "Estimate"]) ^ 2 + 
                        attributes(su.re$varcor$promise)$stddev ^ 2)
      re.weights.cl <- prop.table(exp(re.prom.cl / comb.sd))
      eff.0 <- rep(0, n.promises)
      eff.0[1:group.size] <- group.diff
      sim.output$cluster.sd.ratio[sim] <- sd(data$cluster * su.re$coefficients["cluster", "Estimate"]) / sd(eff.0)
      sim.output$re.cl.sd.ratio[sim] <- attributes(su.re$varcor$promise)$stddev / 0.5
      
    }  else {
      re.weights.cl <- re.weights
      re.prom.cl <- re.prom
      comb.sd <- attributes(su.re.basic$varcor$promise)$stddev
    }
    
    
    
    
    true.eff[sim, ] <- true.values - mean(true.values)
    re.eff[sim, ] <- re.prom - mean(re.prom)
    re.cl.eff[sim, ] <- re.prom.cl - mean(re.prom.cl)
    fe.eff[sim, ] <- fe - mean(fe)
    
    true.wt[sim, ] <- true.weights
    re.wt[sim, ] <- re.weights
    re.cl.wt[sim, ] <- re.weights.cl
    fe.wt[sim, ] <- fe.weights
    
    sim.output$true.sd <- sd(true.values)
    sim.output$fe.sd[sim ] <- sd(fe)
    sim.output$re.sd[sim ] <- attributes(su.re.basic$varcor$promise)$stddev
    sim.output$re.clsd[sim ] <- comb.sd
    
    
    
  }
  
  
  
  
  mean.5.val <- dtf(true = apply(true.eff, 1, function(x) mean(sort(x, decreasing = T)[1:5], na.rm = T)),
                    re = apply(re.eff, 1, function(x) mean(sort(x, decreasing = T)[1:5], na.rm = T)),
                    fe = apply(fe.eff, 1, function(x) mean(sort(x, decreasing = T)[1:5], na.rm = T)),
                    re.cl = apply(re.cl.eff, 1, function(x) mean(sort(x, decreasing = T)[1:5], na.rm = T)))
  
  
  top.5.wt <- dtf(true = apply(true.wt, 1, function(x) sum(sort(x, decreasing = T)[1:5])), 
                  re = apply(re.wt, 1, function(x) sum(sort(x, decreasing = T)[1:5])), 
                  fe = apply(fe.wt, 1, function(x) sum(sort(x, decreasing = T)[1:5])), 
                  re.cl = apply(re.cl.wt, 1, function(x) sum(sort(x, decreasing = T)[1:5])))
  
  sum.diffs <- dtf(re = rowSums(abs(re.wt - true.wt)),
                   fe = rowSums(abs(fe.wt - true.wt)),
                   re.cl = rowSums(abs(re.cl.wt - true.wt)))
  
  
  gini.wt <- dtf(true = apply(true.wt, 1, ineq, type = "Gini"), 
                 re = apply(re.wt, 1, ineq, type = "Gini"), 
                 fe = apply(fe.wt, 1, ineq, type = "Gini"), 
                 re.cl = apply(re.cl.wt, 1, ineq, type = "Gini"))
  
  
  summary.diffs <- dtf(Model = c("RE", "FE", "RE/FMM"),
                       Gap = group.diff,
                       Sims = n.sims,
                       True.sd = mean(sim.output$true.sd),
                       Prop.wt.5 = c(mean(top.5.wt$re - top.5.wt$true, na.rm = T),
                                     mean(top.5.wt$fe - top.5.wt$true, na.rm = T), 
                                     mean(top.5.wt$re.cl - top.5.wt$true, na.rm = T)),
                       mean.wt.5 = c(mean(mean.5.val$re - mean.5.val$true, na.rm = T), 
                                     mean(mean.5.val$fe - mean.5.val$true, na.rm = T), 
                                     mean(mean.5.val$re.cl - mean.5.val$true, na.rm = T)),
                       Gini = c(mean(gini.wt$re - gini.wt$true, na.rm = T),
                                mean(gini.wt$fe - gini.wt$true, na.rm = T), 
                                mean(gini.wt$re.cl - gini.wt$true, na.rm = T)),
                       SD = c(mean(sim.output$re.sd / sim.output$true.sd , na.rm = T),
                              mean(sim.output$fe.sd / sim.output$true.sd , na.rm = T),
                              mean(sim.output$re.clsd / sim.output$true.sd , na.rm = T)), 
                       Abs.Prop.wt.5 = c(mean(abs(top.5.wt$re - top.5.wt$true), na.rm = T),
                                         mean(abs(top.5.wt$fe - top.5.wt$true), na.rm = T), 
                                         mean(abs(top.5.wt$re.cl - top.5.wt$true), na.rm = T)),
                       Abs.mean.wt.5 = c(mean(abs(mean.5.val$re - mean.5.val$true), na.rm = T), 
                                         mean(abs(mean.5.val$fe - mean.5.val$true), na.rm = T), 
                                         mean(abs(mean.5.val$re.cl - mean.5.val$true), na.rm = T)),
                       Abs.Gini = c(mean(abs(gini.wt$re - gini.wt$true), na.rm = T),
                                    mean(abs(gini.wt$fe - gini.wt$true), na.rm = T), 
                                    mean(abs(gini.wt$re.cl - gini.wt$true), na.rm = T)),
                       sum.diff.wt = c(mean(sum.diffs$re, na.rm = T), 
                                       mean(sum.diffs$fe, na.rm = T), 
                                       mean(sum.diffs$re.cl, na.rm = T))
  )
  
  
  # summary.diffs[-1:-3] <- round(summary.diffs[-1:-3], 3)
  
  
  output <- list(weights = list(true = true.weights, re = re.wt, fe = fe.wt, re.cl = re.cl.wt), 
                 effects = list(true= true.eff, re = re.eff, fe = fe.eff, re.cl = re.cl.eff), 
                 summary = summary.diffs, 
                 gini = gini.wt, 
                 top.5.wt = top.5.wt, 
                 mean.5.val = mean.5.val,
                 sim.output = sim.output)
  
  return(output)
}

rmseGroup <- function(x) {
  dtf(re = mean(sqrt(rowMeans((x$effects$re - x$effects$true) ^ 2)), na.rm = T),
      fe = mean(sqrt(rowMeans((x$effects$fe - x$effects$true) ^ 2)), na.rm = T),
      re.cl = mean(sqrt(rowMeans((x$effects$re.cl - x$effects$true) ^ 2)), na.rm = T))
}


rmseGroupWt <- function(x) {
  dtf(re = mean(sqrt(rowMeans((x$weights$re - x$weights$true) ^ 2)), na.rm = T),
      fe = mean(sqrt(rowMeans((x$weights$fe - x$weights$true) ^ 2)), na.rm = T),
      re.cl = mean(sqrt(rowMeans((x$weights$re.cl - x$weights$true) ^ 2)), na.rm = T))
}

set.seed(2139824)
gd5gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 5, group.size = 5, 
                       obs.prom = 100,resp.sd = 1.7)

gd4gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 4, group.size = 5, 
                       obs.prom = 100,resp.sd = 1.7)

gd3gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 3, group.size = 5, 
                       obs.prom = 100,resp.sd = 1.7)


gd2.5gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 2.5, group.size = 5, 
                       obs.prom = 100,resp.sd = 1.7)

gd2gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 2, group.size = 5, 
                       obs.prom = 100,resp.sd = 1.7)

gd1.25gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1.25, group.size = 5, 
                         obs.prom = 100,resp.sd = 1.7)

gd1.5gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1.5, group.size = 5, 
                         obs.prom = 100,resp.sd = 1.7)

gd1gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1, group.size = 5, 
                       obs.prom = 100,resp.sd = 1.7)



gd0.5gs5 <- simDiffsMods(n.sims = 2000, n.promises = 250, 
                         group.diff = 0.5, group.size = 5, 
                         obs.prom = 100,resp.sd = 1.7)

gd5gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 5, group.size = 10, 
                        obs.prom = 100,resp.sd = 1.7)

gd4gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 4, group.size = 10, 
                        obs.prom = 100,resp.sd = 1.7)

gd3gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 3, group.size = 10, 
                        obs.prom = 100,resp.sd = 1.7)

gd2.5gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 2.5, group.size = 10, 
                        obs.prom = 100,resp.sd = 1.7)

gd2gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 2, group.size = 10, 
                        obs.prom = 100,resp.sd = 1.7)

gd1.25gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1.25, group.size = 10, 
                          obs.prom = 100,resp.sd = 1.7)

gd1.5gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1.5, group.size = 10, 
                          obs.prom = 100,resp.sd = 1.7)

gd1gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1, group.size = 10, 
                        obs.prom = 100,resp.sd = 1.7)

gd0.5gs10 <- simDiffsMods(n.sims = 2000, n.promises = 250, 
                          group.diff = 0.5, group.size = 10, 
                          obs.prom = 100,resp.sd = 1.7)




gd5gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 5, group.size = 15, 
                        obs.prom = 100,resp.sd = 1.7)

gd4gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 4, group.size = 15, 
                        obs.prom = 100,resp.sd = 1.7)

gd3gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 3, group.size = 15, 
                        obs.prom = 100,resp.sd = 1.7)

gd2.5gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 2.5, group.size = 15, 
                        obs.prom = 100,resp.sd = 1.7)

gd2gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 2, group.size = 15, 
                        obs.prom = 100,resp.sd = 1.7)
library(mellonMisc)
library(lme4)
library(lfe)
library(mixtools)
library(ineq)
library(reshape2)
set.seed(123898)

gd1.5gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1.5, group.size = 15, 
                          obs.prom = 100,resp.sd = 1.7)


gd1.25gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1.25, group.size = 15, 
                          obs.prom = 100,resp.sd = 1.7)

gd1gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 1, group.size = 15, 
                        obs.prom = 100,resp.sd = 1.7)

# not run?
set.seed(86856)
gd0.5gs15 <- simDiffsMods(n.sims = 2000, n.promises = 250, 
                          group.diff = 0.5, group.size = 15, 
                          obs.prom = 100,resp.sd = 1.7)



gd0 <- simDiffsMods(n.sims = 2000, n.promises = 250, group.diff = 0, group.size = 5, 
                    obs.prom = 100,resp.sd = 1.7)



# load("C:/Dropbox/Manifesto_conjoint/RE_FE_Simulations.Rdata")

library(mellonMisc)
gd3gs5$summary$Gap / gd3gs5$summary["True.sd"]


rmseCalcs <- function(sizes, l.in) {
  rmse.groups <- do.call(rbind, lapply(l.in, rmseGroup))
  rmse.groups <- rmse.groups[, c("fe", "re", "re.cl")]
  colnames(rmse.groups) <- c("FE", "RE", "RE/FMM")
  rmse.groups$Diff <- sizes
  rmse.groups$SD <- sapply(l.in, function(x) x$summary$True.sd[1])
  rmse.groups$gap.sd.ratio  <- rmse.groups$Diff / rmse.groups$SD
  
  rmse.groups <- melt(rmse.groups, id.vars = c("Diff", "gap.sd.ratio", "SD"))
  rmse.groups$Model <- factor(rmse.groups$variable, levels = c("FE", "RE", "RE/FMM"))
  
  rmse.groups.wt <- do.call(rbind, lapply(l.in, rmseGroupWt))
  rmse.groups.wt <- rmse.groups.wt[, c("fe", "re", "re.cl")]
  colnames(rmse.groups.wt) <- c("FE", "RE", "RE/FMM")
  rmse.groups.wt$Diff <- sizes
  rmse.groups.wt$SD <- sapply(l.in, function(x) x$summary$True.sd[1])
  rmse.groups.wt$gap.sd.ratio  <- rmse.groups.wt$Diff / rmse.groups.wt$SD
  
  rmse.groups.wt <- melt(rmse.groups.wt, id.vars = c("Diff", "gap.sd.ratio", "SD"))
  rmse.groups.wt$Model <- factor(rmse.groups.wt$variable, levels = c("FE", "RE", "RE/FMM"))
  
  rmse.wt.plot <- ggplot(data = rmse.groups.wt, aes(x = gap.sd.ratio,
                                                    y = value, group = Model, colour =Model)) + 
    geom_point() + 
    geom_line() +
    scale_y_continuous(limits = c(0, max(rmse.groups.wt$value))) + 
    ylab("RMSE (weights)") + theme_bes()+ xlab("Group gap/SD")+ 
    geom_hline(yintercept = 0, linetype = 2)
  
  
  rmse.plot <- ggplot(data = rmse.groups, aes(x = gap.sd.ratio,
                                              y = value, group = Model, colour =Model)) + 
    geom_point() + 
    geom_line() +
    scale_y_continuous(limits = c(0, max(rmse.groups$value))) + 
    ylab("RMSE (effects)") + theme_bes()+ xlab("Group gap/SD") + 
    geom_hline(yintercept = 0, linetype = 2)
  out <- list(rmse.plot = rmse.plot, rmse.wt.plot = rmse.wt.plot) 
  return(out)
}


library(ggplot2)
comb.size.5 <- list(gd3gs5, gd2gs5, gd1.5gs5, gd1gs5, gd0.5gs5,gd0)
comb.size.10 <- list(gd3gs10, gd2gs10, gd1.5gs10, gd1gs10, gd0.5gs10,gd0)
comb.size.15 <- list(gd3gs15, gd2gs15, gd1.5gs15, gd1gs15, gd0.5gs15,gd0)


rmse.5 <- rmseCalcs(sizes = c(3,2,1.5,1,0.5, 0), l.in = comb.size.5)
rmse.10 <- rmseCalcs(sizes = c(3,2,1.5,1,0.5, 0), l.in = comb.size.10)
rmse.15 <- rmseCalcs(sizes = c(3,2,1.5,1,0.5, 0), l.in = comb.size.15)
library(gridExtra)
rmse.grid<- grid_arrange_shared_legend(rmse.5$rmse.plot + ggtitle("Effects (size=5)"), 
                                       rmse.5$rmse.wt.plot + ggtitle("Weights (size=5)"),
                           rmse.10$rmse.plot + ggtitle("Effects (size=10)"), 
                           rmse.10$rmse.wt.plot + ggtitle("Weights (size=10)"),
                           rmse.15$rmse.plot + ggtitle("Effects (size=15)"), 
                           rmse.15$rmse.wt.plot + ggtitle("Weights (size=15)"), ncol = 2, nrow = 3)
plot(rmse.grid)

saveForPub(rmse.grid, file = "figures/rmse.grid")

rmse.groups <- rbind(rmseGroup(gd3gs5),
                     rmseGroup(gd2gs5),
                     rmseGroup(gd1.5gs5),
                     rmseGroup(gd0.5gs5),
                     rmseGroup(gd1gs5),
                     rmseGroup(gd0))

rmse.groups <- rbind(rmseGroup(gd5gs15),
                     rmseGroup(gd4gs15),
                     rmseGroup(gd3gs15),
                     rmseGroup(gd2gs15),
                     rmseGroup(gd1.5gs15),
                     rmseGroup(gd0.5gs15),
                     rmseGroup(gd1gs15),
                     rmseGroup(gd0))





library(ggplot2)
library(reshape2)
library(tibble)


# this doesn't run properly
levels(rmse.groups$variable)[levels(rmse.groups$variable)=="re"] <- "RE"
levels(rmse.groups$variable)[levels(rmse.groups$variable)=="fe"] <- "FE"
levels(rmse.groups$variable)[levels(rmse.groups$variable)=="re.cl"] <- "RE/FMM"



rmse.groups.wt <- rbind(rmseGroupWt(gd5gs5),
                        rmseGroupWt(gd4gs5),
                        rmseGroupWt(gd3gs5),
                        rmseGroupWt(gd2gs5),
                        rmseGroupWt(gd1.5gs5),
                        rmseGroupWt(gd0.5gs5),
                        rmseGroupWt(gd1gs5),
                        rmseGroupWt(gd0))

rmse.groups.wt$Diff <- c(5,4,3,2,1.5,1,0.5, 0)

# this doesn't run properly
levels(rmse.groups.wt$variable)[levels(rmse.groups.wt$variable)=="re"] <- "RE"
levels(rmse.groups.wt$variable)[levels(rmse.groups.wt$variable)=="fe"] <- "FE"
levels(rmse.groups.wt$variable)[levels(rmse.groups.wt$variable)=="re.cl"] <- "RE/FMM"




comb.sum.5 <- rbind(gd0$summary, gd0.5gs5$summary, gd1gs5$summary, 
                    gd1.5gs5$summary, gd2gs5$summary, gd3gs5$summary, 
                    gd4gs5$summary, gd5gs5$summary)

comb.sum.10 <- rbind(gd0$summary, gd0.5gs10$summary,
                     gd1gs10$summary, gd1.5gs10$summary, 
                     gd2gs10$summary, gd3gs10$summary,
                     gd4gs10$summary, gd5gs10$summary)

comb.sum.15 <- rbind(gd0$summary, gd0.5gs15$summary,
                     gd1gs15$summary, gd1.5gs15$summary, 
                     gd2gs15$summary, gd3gs15$summary,
                     gd4gs15$summary,gd5gs15$summary)

combPlots <- function(comb.sum) {
  
  comb.sum$gap.sd.ratio <- comb.sum$Gap / comb.sum$True.sd
  
  #### Proportion weight top 5 ####
  
  prop.wt5plot <- ggplot(comb.sum, aes(x = gap.sd.ratio, group = Model, colour = Model, y = Prop.wt.5))+ 
    geom_point() + 
    geom_line() +
    # scale_y_continuous(limits = c(0, max(rmse.groups$value))) + 
    theme_bes() + ylab("Mean error in proportion of weight top 5 (bias)") + 
    geom_hline(yintercept = 0, linetype = 2)+ xlab("Group gap/SD")
  
  abs.wt5plot <- ggplot(comb.sum, aes(x = gap.sd.ratio, group = Model, shape= Model, colour = Model, y = Abs.Prop.wt.5))+ 
    geom_point() + 
    geom_line() +
    # scale_y_continuous(limits = c(0, max(rmse.groups$value))) + 
    theme_bes() + ylab("Mean error in proportion of weight top 5 (efficiency)") + 
    geom_hline(yintercept = 0, linetype = 2)+ xlab("Group gap/SD")
  
  #### Standard deviation
  
  sdratioplot <- ggplot(comb.sum, aes(x = gap.sd.ratio, group = Model, colour = Model, shape= Model, y = SD))+ 
    geom_point() + 
    geom_line() +
    # scale_y_continuous(limits = c(0, max(rmse.groups$value))) + 
    ylab("Ratio of estimated effect SD to true effect SD") + theme_bes()  + 
    geom_hline(yintercept = 1, linetype = 2)+ xlab("Group gap/SD")
  
  
  
  #### Gini ####
  gini.plot <- ggplot(comb.sum, aes(x = gap.sd.ratio, group = Model, shape= Model, colour = Model, y = Gini))+ 
    geom_point() + 
    geom_line() +
    # scale_y_continuous(limits = c(0, max(rmse.groups$value))) + 
    theme_bes() + ylab("Mean error Gini") + geom_hline(yintercept = 0, linetype = 2)+ xlab("Group gap/SD")
  
  
  absgini.plot <- ggplot(comb.sum, aes(x = gap.sd.ratio, group = Model, shape= Model, colour = Model, y = Abs.Gini))+ 
    geom_point() + 
    geom_line() +
    # scale_y_continuous(limits = c(0, max(rmse.groups$value))) + 
    ylab("Mean absolute error Gini") + theme_bes()  + 
    geom_hline(yintercept = 0, linetype = 2) + xlab("Group gap/SD")
  
  #### Difference weights
  
  diff.wt.plot <- ggplot(comb.sum, aes(x = gap.sd.ratio, group = Model, colour = Model, y = sum.diff.wt))+ 
    geom_point() + 
    geom_line() +
    # scale_y_continuous(limits = c(0, max(rmse.groups$value))) + 
    ylab("Sum of differences in weights (max value 2)") + theme_bes()  + 
    geom_hline(yintercept = 0, linetype = 2)+ xlab("Group gap/SD")
  
  
  output <- list(prop.wt5plot = prop.wt5plot, 
                 abs.wt5plot = abs.wt5plot,
                 sdratioplot = sdratioplot,
                 gini.plot = gini.plot,
                 absgini.plot = absgini.plot,
                 diff.wt.plot = diff.wt.plot)
  return(output)
}

comb.plot.5 <- combPlots(comb.sum.5)
comb.plot.10 <- combPlots(comb.sum.10)
comb.plot.15 <- combPlots(comb.sum.15)


library(grid)
library(gridExtra)


gini.error.plot <- grid_arrange_shared_legend(comb.plot.5$gini.plot + ggtitle("Bias (Size=5)")+ ylab("Mean error Gini"),  
                                comb.plot.5$absgini.plot + ggtitle("Efficiency (Size=5)") + ylab("Mean absolute error Gini"), 
                                comb.plot.10$gini.plot  + ggtitle("Bias (Size=10)")+ ylab("Mean error Gini"),
                                comb.plot.10$absgini.plot  + ggtitle("Efficiency (Size=10)") + ylab("Mean absolute error Gini"), 
                                comb.plot.15$gini.plot + ggtitle("Bias (Size=15)") + ylab("Mean error Gini"), 
                                comb.plot.15$absgini.plot + ggtitle("Efficiency (Size=15)") + ylab("Mean absolute error Gini"),
                                ncol = 2, nrow = 3)

plot(gini.error.plot)

saveForPub(gini.error.plot, file = "figures/gini.error.plot")

top5error.plot <- grid_arrange_shared_legend(comb.plot.5$prop.wt5plot + ggtitle("Bias (Size=5)")+ ylab("Mean error top 5 weight"),  
                                comb.plot.5$abs.wt5plot + ggtitle("Efficiency (Size=5)") + ylab("Mean absolute error top 5 weight"), 
                                comb.plot.10$prop.wt5plot  + ggtitle("Bias (Size=10)")+ ylab("Mean error top 5 weight"),
                                comb.plot.10$abs.wt5plot  + ggtitle("Efficiency (Size=10)") + ylab("Mean absolute error top 5 weight"), 
                                comb.plot.15$prop.wt5plot + ggtitle("Bias (Size=15)") + ylab("Mean error top 5 weight"), 
                                comb.plot.15$abs.wt5plot + ggtitle("Efficiency (Size=15)") + ylab("Mean absolute error top 5 weight"),
                                ncol = 2, nrow = 3)

saveForPub(top5error.plot, file = "figures/top5error.plot")


true.v.observed.sd <- grid_arrange_shared_legend(comb.plot.5$sdratioplot + ggtitle("Size=5)")+ ylab("Observed SD/True SD"),  
                                comb.plot.10$sdratioplot  + ggtitle("Size=10")+ ylab("Observed SD/True SD"),
                                comb.plot.15$sdratioplot + ggtitle("Size=15") + ylab("Observed SD/True SD"), 
                                ncol = 1, nrow = 3)


plot(true.v.observed.sd)
saveForPub(true.v.observed.sd, file = "figures/true.v.observed.sd")

rmse.groups$Diff

comb.sum.5

comb.plot.10$abs.wt5plot

comb.plot.15$prop.wt5plot
comb.plot.15$abs.wt5plot



comb.plot.10$sdratioplot


comb.plot.5$diff.wt.plot


comb.plot.10$diff.wt.plot

comb.plot.15$diff.wt.plot


save.image("data/intermediate/RE_FE_Simulations.Rdata")
